import numpy as np
import matplotlib.pyplot as plt

from activefluids.Module.defects import *

data_path = r'D:\Scientific_data\Experimental_data_defects'
plot_path = r'D:\Scientific_data\Plots'

Check for experimental example

data_path = r'D:\Scientific_data\Experimental_data_defects'
plot_path = r'D:\Scientific_data\Plots'

# file = r'\4.jpg'
file = r'\1-40-100\1-0001.tif'

img = plt.imread(data_path + file)

# from skimage.color import rgb2gray
# img = rgb2gray(img)

shift = 1/2 # Due to the choice of the 2x2 grid (neigbors in the northeastern corner), we move the detected defect to the middle of the cell

pix_x = img.shape[1]
pix_y = img.shape[0]

x = np.arange(0,pix_x)
y = np.arange(0,pix_y)

xx, yy = np.meshgrid(x, y)

Giomi

ori, coh, E = orientation_analysis(img, 15)
k = compute_topological_charges(ori, int_area='cell', origin='upper')
defects = localize_defects(k, x_grid=xx, y_grid=yy)
compute_defect_orientations(ori, defects)

plushalf = defects[defects['charge']==.5]
minushalf = defects[defects['charge']==-.5]

fig, ax  = plt.subplots(figsize=(16,16))
s = 15
ax.imshow(img,cmap='gray')
ax.quiver(xx[::s,::s], yy[::s,::s], np.cos(ori)[::s,::s], np.sin(ori)[::s,::s], headaxislength=0, headwidth=0, headlength=0, color='y', scale=60, pivot='mid')

ax.plot(plushalf['x'], plushalf['y'],'ro',markersize=15,label=r'+1/2 defect')
ax.quiver(plushalf['x'], plushalf['y'], np.cos(plushalf['ang1']), np.sin(plushalf['ang1']), headaxislength=0, headwidth=0, headlength=0, color='r', scale=50)

for i in range(3):
    ax.quiver(minushalf['x'], minushalf['y'], np.cos(minushalf['ang'+str(i+1)]), np.sin(minushalf['ang'+str(i+1)]), headaxislength=0, headwidth=0, headlength=0, color='b', scale=50)

ax.set_xlabel('x (in pixels)')
ax.set_ylabel('y (in pixels)')

Text(0, 0.5, 'y (in pixels)')

Simple

The image (and all subsequent computations) have to be flipped in order to work properly. Maybe this can be streamlined/avoided by some additional keywords??

ori, coh, E = orientation_analysis(img, 15)
k = compute_topological_charges(-ori, int_area='cell', origin='lower')
defects = localize_defects(k, x_grid=xx, y_grid=yy)
compute_defect_orientations(-ori, defects, method='interpolation', x_grid=x, y_grid=y, interpolation_radius=5,  min_sep=1)

plushalf = defects[defects['charge']==.5]
minushalf = defects[defects['charge']==-.5]

fig, ax  = plt.subplots(figsize=(16,16))
s = 15
ax.imshow(img,cmap='gray',origin='lower')
ax.quiver(xx[::s,::s], yy[::s,::s], np.cos(ori)[::s,::s], -np.sin(ori)[::s,::s], headaxislength=0, headwidth=0, headlength=0, color='y', scale=60, pivot='mid')

ax.plot(plushalf['x'], plushalf['y'],'ro',markersize=15,label=r'+1/2 defect')
ax.quiver(plushalf['x'], plushalf['y'], np.cos(plushalf['ang1']), np.sin(plushalf['ang1']), headaxislength=0, headwidth=0, headlength=0, color='r', scale=50)

for i in range(3):
    ax.quiver(minushalf['x'], minushalf['y'], np.cos(minushalf['ang'+str(i+1)]), np.sin(minushalf['ang'+str(i+1)]), headaxislength=0, headwidth=0, headlength=0, color='b', scale=50)

ax.set_xlabel('x (in pixels)')
ax.set_ylabel('y (in pixels)')

Text(0, 0.5, 'y (in pixels)')

# fig.savefig(r'D:\Scientific_data\Plots\defect_orientation_asymetric.png', bbox_inches = 'tight', dpi=100)

diff_1 = modulo(minushalf['ang1']-minushalf['ang2'],type='polar')
diff_2 = modulo(minushalf['ang1']-minushalf['ang3'],type='polar')
diff_3 = modulo(minushalf['ang2']-minushalf['ang3'],type='polar')
diff = np.hstack((diff_1,diff_2,diff_3))

np.rad2deg(diff)

array([-115.2,  111.6,  -93.6, -118.8,   97.2,  100.8, -111.6,  104.4,
       -115.2, -118.8,  108. ,  -97.2, -122.4, -129.6, -115.2,  108. ,
        104.4, -118.8, -111.6,  122.4,  111.6,  115.2, -118.8,  108. ,
        129.6,  115.2,   97.2,  111.6, -144. ,  136.8,  108. , -126. ,
       -122.4,  126. , -122.4,  122.4,  122.4, -122.4,  129.6,  122.4,
        118.8,  129.6, -118.8, -136.8,  122.4,  129.6, -126. , -115.2,
       -122.4,  104.4, -133.2, -108. , -115.2, -118.8, -133.2,  104.4,
       -129.6, -133.2,  136.8,  136.8, -122.4,  133.2, -122.4, -118.8,
        129.6, -133.2, -115.2, -111.6, -115.2,  133.2,  118.8, -118.8,
       -118.8,  111.6,  133.2,  122.4, -136.8,  118.8,  122.4,  129.6,
        144. ])

fig, ax = plt.subplots(figsize=(5,5))
ax.hist(np.abs(np.rad2deg(diff)),bins='auto',density=True)
ax.set_xlabel('Degree difference between -1/2 orientations')
ax.set_ylabel('PDF')

Text(0, 0.5, 'PDF')

print(np.argmin(diff_1))
print(np.argmin(diff_2))
print(np.argmin(diff_3))

13
1
22

print(np.rad2deg(diff_1[1]))
print(np.rad2deg(diff_2[1]))
print(np.rad2deg(diff_3[1]))

111.59999999999998
-144.00000000000006
104.39999999999996

# fig, ax  = plt.subplots(figsize=(16,16))
# s = 15
# ax.imshow(img,cmap='gray',origin='lower')
# ax.quiver(xx[::s,::s], yy[::s,::s], np.cos(ori)[::s,::s], np.sin(ori)[::s,::s], headaxislength=0, headwidth=0, headlength=0, color='y', scale=60, pivot='mid')

# ax.plot(plushalf['x'], plushalf['y'],'ro',markersize=15,label=r'+1/2 defect')
# ax.quiver(plushalf['x'], plushalf['y'], np.cos(plushalf['ang1']), np.sin(plushalf['ang1']), headaxislength=0, headwidth=0, headlength=0, color='r', scale=50)

# for i in range(3):
#     ax.quiver(minushalf['x'], minushalf['y'], np.cos(minushalf['ang'+str(i+1)]), np.sin(minushalf['ang'+str(i+1)]), headaxislength=0, headwidth=0, headlength=0, color='b', scale=50)

# ax.set_xlabel('x (in pixels)')
# ax.set_ylabel('y (in pixels)')

Check consistency of quiver/imshow

data_path = r'D:\Scientific_data\Experimental_data_defects'
plot_path = r'D:\Scientific_data\Plots'

file = r'\4.jpg'
# file = r'\1-40-100\1-0001.tif'

img = plt.imread(data_path + file)

from skimage.color import rgb2gray
img = rgb2gray(img)

shift = 1/2 # Due to the choice of the 2x2 grid (neigbors in the northeastern corner), we move the detected defect to the middle of the cell

pix_x = img.shape[1]
pix_y = img.shape[0]

x = np.arange(0,pix_x)
y = np.arange(0,pix_y)

xx, yy = np.meshgrid(x, y)

ori, coh, E = orientation_analysis(img, 15)
k = compute_topological_charges(ori, int_area='cell', origin='upper')
defects = localize_defects(k, x_grid=xx, y_grid=yy)
compute_defect_orientations(ori, defects)

fig, axs  = plt.subplots(1,2,figsize=(24,12))
s = 15
axs[0].imshow(img,cmap='gray')
axs[0].quiver(xx[::s,::s], yy[::s,::s], np.cos(ori)[::s,::s], np.sin(ori)[::s,::s], headaxislength=0, headwidth=0, headlength=0, color='r', scale=60, pivot='mid')
axs[1].imshow(img,cmap='gray',origin='lower')
axs[1].quiver(xx[::s,::s], yy[::s,::s], np.cos(-ori)[::s,::s], np.sin(-ori)[::s,::s], headaxislength=0, headwidth=0, headlength=0, color='r', scale=60, pivot='mid')

<matplotlib.quiver.Quiver at 0x29f0d6ced40>

point_x = 600
point_y = 400
ang = 3*np.pi/4

fig, axs  = plt.subplots(1,2,figsize=(24,12))

s = 15
axs[0].imshow(img,cmap='gray')
axs[0].quiver(point_x, point_y, np.cos(ang), np.sin(ang), color='r', scale=6)
axs[1].imshow(img,cmap='gray',origin='lower')
axs[1].quiver(point_x, point_y, np.cos(-ang), np.sin(-ang), color='r', scale=6)

<matplotlib.quiver.Quiver at 0x29f0110ffd0>

arr = np.zeros((3,3))
arr[2,1] = 1
print(arr)

[[0. 0. 0.]
 [0. 0. 0.]
 [0. 1. 0.]]

Analytical examples

Symmetric defects as minimizers of elastic energy (equal Frank constants)

fig, axs = plt.subplots(3,2,figsize=(12,18),sharex=True,sharey=True)

axs[0,0].set_title('Direct method')
axs[0,1].set_title('Giomi')

charges = [1/2,-1/2,-1/2]
theta_0 = [np.pi/3,np.pi/3,0]

x = np.arange(-5,5,.5)
y = np.arange(-5,5,.5)
xx, yy = np.meshgrid(x,y)

phi = np.arctan2(yy,xx)

for i in range(len(charges)):
    theta = charges[i]*phi + theta_0[i]

    k = compute_topological_charges(theta, int_area='square', width=1, origin='lower')
    defects = localize_defects(k, x_grid=xx, y_grid=yy)
    defects['x'] = 0
    defects['y'] = 0
    for j in range(2):
        axs[i,j].quiver(xx, yy, np.cos(theta), np.sin(theta), headaxislength=0, headwidth=0, headlength=0, pivot='mid')

        if j==0:
            compute_defect_orientations(theta, defects, method='interpolation', x_grid=xx, y_grid=yy, min_sep=1, interpolation_radius=1)
        else:
            compute_defect_orientations(theta, defects)

        plushalf = defects[defects['charge']==.5]
        minushalf = defects[defects['charge']==-.5]

        axs[i,j].plot(plushalf['x'], plushalf['y'],'ro',markersize=15,label=r'+1/2 defect')
        axs[i,j].plot(plushalf['x'], plushalf['y'],'ro',markersize=15,label=r'+1/2 defect')
        axs[i,j].quiver(plushalf['x'], plushalf['y'], np.cos(plushalf['ang1']), np.sin(plushalf['ang1']), headaxislength=0, headwidth=0, headlength=0, color='r', scale=5)

        if 'ang2' in defects:
            for k in range(3):
                axs[i,j].quiver(minushalf['x'], minushalf['y'], np.cos(minushalf['ang'+str(k+1)]), np.sin(minushalf['ang'+str(k+1)]), headaxislength=0, headwidth=0, headlength=0, color='b', scale=6)
                axs[i,j].quiver(minushalf['x'], minushalf['y'], np.cos(minushalf['ang'+str(k+1)]), np.sin(minushalf['ang'+str(k+1)]), headaxislength=0, headwidth=0, headlength=0, color='b', scale=2)

C:\Users\Vasco\anaconda3\envs\default-env\lib\site-packages\scipy\interpolate\_fitpack_impl.py:977: RuntimeWarning: No more knots can be added because the number of B-spline
coefficients already exceeds the number of data points m.
Probable causes: either s or m too small. (fp>s)
	kx,ky=1,1 nx,ny=24,21 m=400 fp=0.000004 s=0.000000
  warnings.warn(RuntimeWarning(_iermess2[ierm][0] + _mess))

Build asymmetric defect by hand

N = 100
R = 1

theta = np.linspace(0,2*np.pi,N,endpoint=False)
x_c = R*np.cos(theta)
y_c = R*np.sin(theta)
delta = .25
L = np.pi-.2

x_grid = np.arange(-L,L,delta)
y_grid = np.arange(-L,L,delta)
xx, yy = np.meshgrid(x_grid,y_grid)

u_lattice = np.sin(yy)
v_lattice = np.sin(xx)
u = np.copy(u_lattice)
v = np.copy(v_lattice)
v[yy<=0] = 1
u[yy<=0] = 0

phi = get_orientation_angle(u, v)
k = compute_topological_charges(phi, int_area='square', width=1, origin='lower')

C:\Users\Vasco\Code\Notebooks\activefluids\Module\defects.py:44: RuntimeWarning: divide by zero encountered in true_divide
  phi = np.arctan(v/u) # TO DO: catch zero exception!!

fig, ax = plt.subplots(figsize=(9,8))
c = ax.pcolormesh(xx,yy,k,shading='auto')
ax.quiver(xx, yy, u, v, headaxislength=0, headwidth=0, headlength=0, pivot='mid')
plt.colorbar(c, label='topological charge')

<matplotlib.colorbar.Colorbar at 0x29f0bfd75b0>

defects = localize_defects(k, x_grid=xx, y_grid=yy)
defects['x'] = 0
defects['y'] = 0
compute_defect_orientations(phi, defects, method='interpolation', x_grid=xx, y_grid=yy, min_sep=1, interpolation_radius=R)
# compute_defect_orientations(phi, defects, origin='lower')
plushalf = defects[defects['charge']==.5]
minushalf = defects[defects['charge']==-.5]

C:\Users\Vasco\anaconda3\envs\default-env\lib\site-packages\scipy\interpolate\_fitpack_impl.py:977: RuntimeWarning: A theoretically impossible result when finding a smoothing spline
with fp = s. Probable causes: s too small or badly chosen eps.
(abs(fp-s)/s>0.001)
	kx,ky=1,1 nx,ny=26,22 m=576 fp=0.000000 s=0.000000
  warnings.warn(RuntimeWarning(_iermess2[ierm][0] + _mess))

defects

	charge	x	y	x_ind	y_ind	ang1	ang2	ang3
0	-0.5	0	0	11	12	-1.570796	0.816814	2.324779

x0, y0 = 0,0
fig, axs = plt.subplots(1,2,figsize=(16,8))
axs[0].quiver(xx, yy, u, v, headaxislength=0, headwidth=0, headlength=0, pivot='mid')

axs[0].plot(plushalf['x'], plushalf['y'],'ro',markersize=15,label=r'+1/2 defect')
axs[0].plot(plushalf['x'], plushalf['y'],'ro',markersize=15,label=r'+1/2 defect')
axs[0].quiver(plushalf['x'], plushalf['y'], np.cos(plushalf['ang1']), np.sin(plushalf['ang1']), headaxislength=0, headwidth=0, headlength=0, color='r', scale=50)

for i in range(3):
    axs[0].quiver(minushalf['x'], minushalf['y'], np.cos(minushalf['ang'+str(i+1)]), np.sin(minushalf['ang'+str(i+1)]), headaxislength=0, headwidth=0, headlength=0, color='b', scale=6)
    axs[1].quiver(minushalf['x'], minushalf['y'], np.cos(minushalf['ang'+str(i+1)]), np.sin(minushalf['ang'+str(i+1)]), headaxislength=0, headwidth=0, headlength=0, color='b', scale=2)

axs[0].plot(x0+x_c,y0+y_c,'b*')
axs[0].plot(x0,y0,'rx')
axs[1].plot(x0+x_c,y0+y_c,'b*')

[<matplotlib.lines.Line2D at 0x29f0c73e830>]

phi = get_orientation_angle(u, -v)
k = compute_topological_charges(phi, int_area='square', width=1, origin='lower')

C:\Users\Vasco\Code\Notebooks\activefluids\Module\defects.py:44: RuntimeWarning: divide by zero encountered in true_divide
  phi = np.arctan(v/u) # TO DO: catch zero exception!!

fig, ax = plt.subplots(figsize=(9,8))
c = ax.pcolormesh(xx,yy,k,shading='auto')
ax.quiver(xx, yy, u, -v, headaxislength=0, headwidth=0, headlength=0, pivot='mid')
plt.colorbar(c, label='topological charge')

<matplotlib.colorbar.Colorbar at 0x29f0be62020>

defects = localize_defects(k, x_grid=xx, y_grid=yy)
defects['x'] = 0
defects['y'] = 0
compute_defect_orientations(phi, defects, method='interpolation', x_grid=xx, y_grid=yy, min_sep=1, interpolation_radius=R)
# compute_defect_orientations(phi, defects, origin='lower')
plushalf = defects[defects['charge']==.5]
minushalf = defects[defects['charge']==-.5]

C:\Users\Vasco\anaconda3\envs\default-env\lib\site-packages\scipy\interpolate\_fitpack_impl.py:977: RuntimeWarning: A theoretically impossible result when finding a smoothing spline
with fp = s. Probable causes: s too small or badly chosen eps.
(abs(fp-s)/s>0.001)
	kx,ky=1,1 nx,ny=26,22 m=576 fp=0.000000 s=0.000000
  warnings.warn(RuntimeWarning(_iermess2[ierm][0] + _mess))

defects

	charge	x	y	x_ind	y_ind	ang1
0	0.5	0	0	11	12	-1.570796

x0, y0 = 0,0
fig, axs = plt.subplots(1,2,figsize=(16,8))
axs[0].quiver(xx, yy, u, -v, headaxislength=0, headwidth=0, headlength=0, pivot='mid')

axs[0].plot(plushalf['x'], plushalf['y'],'ro',markersize=15,label=r'+1/2 defect')
axs[1].plot(plushalf['x'], plushalf['y'],'ro',markersize=15,label=r'+1/2 defect')
axs[0].quiver(plushalf['x'], plushalf['y'], np.cos(plushalf['ang1']), np.sin(plushalf['ang1']), headaxislength=0, headwidth=0, headlength=0, color='r', scale=5)
axs[1].quiver(plushalf['x'], plushalf['y'], np.cos(plushalf['ang1']), np.sin(plushalf['ang1']), headaxislength=0, headwidth=0, headlength=0, color='r', scale=3)

# for i in range(3):
#     axs[0].quiver(minushalf['x'], minushalf['y'], np.cos(minushalf['ang'+str(i+1)]), np.sin(minushalf['ang'+str(i+1)]), headaxislength=0, headwidth=0, headlength=0, color='b', scale=6)
#     axs[1].quiver(minushalf['x'], minushalf['y'], np.cos(minushalf['ang'+str(i+1)]), np.sin(minushalf['ang'+str(i+1)]), headaxislength=0, headwidth=0, headlength=0, color='b', scale=2)

axs[0].plot(x0+x_c,y0+y_c,'b*')
axs[0].plot(x0,y0,'rx')
axs[1].plot(x0+x_c,y0+y_c,'b*')

[<matplotlib.lines.Line2D at 0x29f0d519f00>]

N = 100
R = 1

theta = np.linspace(0,2*np.pi,N,endpoint=False)
x_c = R*np.cos(theta)
y_c = R*np.sin(theta)

delta = .25
L = np.pi-.2

x_grid = np.arange(-L,L,delta)
y_grid = np.arange(-L,L,delta)
xx, yy = np.meshgrid(x_grid,y_grid)

u_lattice = np.sin(yy)
v_lattice = np.sin(xx)
u = np.copy(u_lattice)
v = np.copy(v_lattice)
v[yy<=xx] = 1
u[yy<=xx] = 1

phi = get_orientation_angle(u, v)
k = compute_topological_charges(phi, int_area='square', width=1, origin='lower')

fig, ax = plt.subplots(figsize=(9,8))
c = ax.pcolormesh(xx,yy,k,shading='auto')
ax.quiver(xx, yy, u, v, headaxislength=0, headwidth=0, headlength=0, pivot='mid')
plt.colorbar(c, label='topological charge')

<matplotlib.colorbar.Colorbar at 0x29f10a27dc0>

defects = localize_defects(k, x_grid=xx, y_grid=yy)
defects['x'] = 0
defects['y'] = 0
compute_defect_orientations(phi, defects, method='interpolation', x_grid=xx, y_grid=yy, min_sep=1, interpolation_radius=R)
# compute_defect_orientations(phi, defects, origin='lower')
plushalf = defects[defects['charge']==.5]
minushalf = defects[defects['charge']==-.5]

C:\Users\Vasco\anaconda3\envs\default-env\lib\site-packages\scipy\interpolate\_fitpack_impl.py:977: RuntimeWarning: No more knots can be added because the number of B-spline
coefficients already exceeds the number of data points m.
Probable causes: either s or m too small. (fp>s)
	kx,ky=1,1 nx,ny=26,27 m=576 fp=0.002241 s=0.000000
  warnings.warn(RuntimeWarning(_iermess2[ierm][0] + _mess))
C:\Users\Vasco\anaconda3\envs\default-env\lib\site-packages\scipy\interpolate\_fitpack_impl.py:977: RuntimeWarning: No more knots can be added because the number of B-spline
coefficients already exceeds the number of data points m.
Probable causes: either s or m too small. (fp>s)
	kx,ky=1,1 nx,ny=26,27 m=576 fp=0.000028 s=0.000000
  warnings.warn(RuntimeWarning(_iermess2[ierm][0] + _mess))

x0, y0 = 0,0
fig, axs = plt.subplots(1,2,figsize=(16,8))
axs[0].quiver(xx, yy, u, v, headaxislength=0, headwidth=0, headlength=0, pivot='mid')

axs[0].plot(plushalf['x'], plushalf['y'],'ro',markersize=15,label=r'+1/2 defect')
axs[0].plot(plushalf['x'], plushalf['y'],'ro',markersize=15,label=r'+1/2 defect')
axs[0].quiver(plushalf['x'], plushalf['y'], np.cos(plushalf['ang1']), np.sin(plushalf['ang1']), headaxislength=0, headwidth=0, headlength=0, color='r', scale=50)

for i in range(3):
    axs[0].quiver(minushalf['x'], minushalf['y'], np.cos(minushalf['ang'+str(i+1)]), np.sin(minushalf['ang'+str(i+1)]), headaxislength=0, headwidth=0, headlength=0, color='b', scale=6)
    axs[1].quiver(minushalf['x'], minushalf['y'], np.cos(minushalf['ang'+str(i+1)]), np.sin(minushalf['ang'+str(i+1)]), headaxislength=0, headwidth=0, headlength=0, color='b', scale=2)

axs[0].plot(x0+x_c,y0+y_c,'b*')
axs[0].plot(x0,y0,'rx')
axs[1].plot(x0+x_c,y0+y_c,'b*')

[<matplotlib.lines.Line2D at 0x29f11211810>]